"""FIle to read a gmsh msh file using FiPy's gmshImport"""

import numpy
import numpy.ma as MA
import tempfile

from grid import MeshType
from matvec import Vertex, Vector

from pylab import find

##########################################################################
# The following code is borrowed from FiPy
##########################################################################
class GmshFile:
    def __init__(self, filename):
        self.filename = filename
        
        self.formatWritten = False
        
        # open the .msh file
        if (isinstance(self.filename, file) 
            or (hasattr(self.filename, "file") 
                and isinstance(self.filename.file, file))):
            self.fileobj = self.filename
            self.filename = self.fileobj.name
        else:
            self.fileobj = open(name=self.filename, mode='r')

    def _getElementType(self, vertices, dimensions):
        if (vertices == 3 and dimensions == 2):
            return 2 ## triangle
        elif (vertices == 4 and dimensions == 2):
            return 3 ## quadrangle
        elif (vertices == 4 and dimensions == 3):
            return 4 ## tetrahedron
        elif (vertices == 8 and dimensions == 3):
            return 5 ## hexahedron
        elif (vertices == 6 and dimensions == 3):
            return 6 ## prism
        elif (vertices == 5 and dimensions == 3):
            return 7 ## pyramid
        else:
            raise MeshExportError, "Element type unsupported by Gmsh"
    
    def _orderVertices(self, vertexCoords, vertices):
        coordinates = numpy.take(vertexCoords, vertices, axis=1)
        centroid = numpy.add.reduce(coordinates, axis=1) / coordinates.shape[1]
        coordinates = coordinates - centroid[..., numpy.newaxis]

        # to prevent div by zero
        coordinates = numpy.where(coordinates == 0, 1.e-10, coordinates)

        # angles go from -pi / 2 to 3*pi / 2
        angles = numpy.arctan(coordinates[1] / coordinates[0]) \
                + numpy.where(coordinates[0] < 0, numpy.pi, 0) 
        sortorder = numpy.argsort(angles)
        return numpy.take(vertices, sortorder)
        
    def write(self, obj, time=0.0, timeindex=0):
        pass
        
    def close(self):
        self.fileobj.close()

class MSHFile(GmshFile):
    """
    Class responsible for parsing a Gmsh file and then readying
    its contents for use by a `Mesh` constructor. 
    
    Can handle a partitioned mesh based on `parallel.Nproc`. If partitioning,
    the msh file must either be previously partitioned with the number of
    partitions matching `Nproc`, or the mesh must be specified with a geo file
    or multiline string.

    Does not support gmsh versions < 2. If partitioning, gmsh
    version must be >= 2.5.

    TODO: Refactor face extraction functions.
    """
    def __init__(self, filename, dimensions, coordDimensions=None):
        """
        :Parameters:
          - `filename`: a string indicating gmsh output file
          - `dimensions`: an integer indicating dimension of mesh
          - `coordDimensions`: an integer indicating dimension of shapes

            The file will be created if it doesn't exist when opened
            for writing; it will be truncated when opened for writing.
            Add a 'b' to the mode for binary files.
        """
        self.dimensions = dimensions
        self.coordDimensions = coordDimensions
        
        self.mesh = None
        self.meshWritten = False

        GmshFile.__init__(self, filename=filename)
        
    def _getMetaData(self):
        """
        Extracts gmshVersion, file-type, and data-size in that
        order.
        """
        self._seekForHeader("MeshFormat")
        metaData = self.fileobj.readline().split()
        self.fileobj.seek(0)
        return [float(x) for x in metaData]

    def _isolateData(self, title):
        """
        Gets all data between $[title] and $End[title], writes
        it out to its own file.
        """
        newF = tempfile.TemporaryFile()
        self._seekForHeader(title)
        
        # extract the actual data within section
        while True:
            line = self.fileobj.readline()
            if ("$End%s" % title) not in line: 
                newF.write(line) 
            else: break

        self.fileobj.seek(0) 
        newF.seek(0) # restore file positions
        return newF

    def _seekForHeader(self, title):
        """
        Iterate through a file until we end up at the section header
        for `title`. Function has obvious side-effects on `self.fileobj`.
        """
        while True:
            line = self.fileobj.readline()
            if len(line) == 0:
                raise EOFError("No `%s' header found!" % title)
                break
            elif (("$%s" % title) not in line): 
                continue
            else: 
                break # found header

    def _deriveCellsAndFaces(self, cellsToVertIDs, shapeTypes, numCells):
        """
        Uses element information obtained from `_parseElementFile` to deliver
        `facesToVertices` and `cellsToFaces`.
        """

        #def formatForFiPy(arr): return arr.swapaxes(0,1)[::-1]
        def formatForFiPy(arr): return arr.swapaxes(0,1)

        allShapes  = numpy.unique(shapeTypes).tolist()
        maxFaces   = max([self.numFacesPerCell[x] for x in allShapes])

        # `cellsToFaces` must be padded with -1; see mesh.py
        currNumFaces = 0
        cellsToFaces = numpy.ones((numCells, maxFaces), 'l') * -1
        facesDict    = {}
        uniqueFaces  = []

        # we now build `cellsToFaces` and `uniqueFaces`,
        # the latter will result in `facesToVertices`.
        for cellIdx in range(numCells):
            shapeType = shapeTypes[cellIdx]
            cell = cellsToVertIDs[cellIdx]

            if shapeType in [5, 12, 17]: # hexahedron
                faces = self._extractOrderedFaces(cell=cell, 
                                                  faceOrderings=[[0, 1, 2, 3], # ordering of vertices gleaned from
                                                                 [4, 5, 6, 7], # a one-cube Grid3D example
                                                                 [0, 1, 5, 4],
                                                                 [3, 2, 6, 7],
                                                                 [0, 3, 7, 4],
                                                                 [1, 2, 6, 5]])
            elif shapeType in [6, 13, 18]: # prism
                faces = self._extractOrderedFaces(cell=cell, 
                                                  faceOrderings=[[0, 1, 2],
                                                                 [5, 4, 3],
                                                                 [3, 4, 1, 0],
                                                                 [4, 5, 2, 1],
                                                                 [5, 3, 0, 2]])
            elif shapeType in [7, 14, 19]: # pyramid
                faces = self._extractOrderedFaces(cell=cell, 
                                                  faceOrderings=[[0, 1, 2, 3],
                                                                 [0, 1, 4],
                                                                 [1, 2, 4],
                                                                 [2, 3, 4],
                                                                 [3, 0, 4]])
            else:
                if shapeType in [2, 9, 20, 21, 22, 23, 24, 25]:
                    faceLength = 2 # triangle
                elif shapeType in [3, 10, 16]:
                    faceLength = 2 # quadrangle
                elif shapeType in [4, 11, 29, 30, 31]:
                    faceLength = 3 # tetrahedron
                
                faces = self._extractRegularFaces(cell=cell, 
                                                  faceLength=faceLength, 
                                                  facesPerCell=self.numFacesPerCell[shapeType])

            for faceIdx in range(len(faces)):
                # NB: currFace is sorted for the key to spot duplicates
                currFace = faces[faceIdx]
                keyStr   = ' '.join([str(x) for x in sorted(currFace)])

                if facesDict.has_key(keyStr):
                    cellsToFaces[cellIdx][faceIdx] = facesDict[keyStr]
                else: # new face
                    facesDict[keyStr] = currNumFaces
                    cellsToFaces[cellIdx][faceIdx] = currNumFaces
                    uniqueFaces.append(currFace)
                    currNumFaces += 1
                    
        # pad short faces with -1
        maxFaceLen = max([len(f) for f in uniqueFaces])
        uniqueFaces = [[-1] * (maxFaceLen - len(f)) + f for f in uniqueFaces]
               
        facesToVertices = numpy.array(uniqueFaces, dtype=int)

        return formatForFiPy(facesToVertices), formatForFiPy(cellsToFaces), facesDict
        #return facesToVertices, cellsToFaces, facesDict

    def _translateNodesToVertices(self, entitiesNodes, vertexMap):
        """Translates entitiesNodes from Gmsh node IDs to `vertexCoords` indices.
        """
        entitiesVertices = []

        for entity in entitiesNodes:
            try:
                vertIndices = vertexMap[numpy.array(entity)]
            except IndexError:
                vertIndices = numpy.ones((len(entity),), 'l') * -1
            entitiesVertices.append(vertIndices)

        return entitiesVertices

    def _extractRegularFaces(self, cell, faceLength, facesPerCell):
        """Return faces for a regular poly(gon|hedron)
        
        Given `cell`, defined by node IDs, returns an array of
        `facesPerCell` faces of length `faceLength` in terms of nodes.
        """
        faces = []
        for i in range(facesPerCell):
            aFace = []
            for j in range(faceLength):
                aVertex = (i + j) % len(cell) # we may wrap
                aFace.append(int(cell[aVertex]))
            faces.append(aFace)
        return faces

    def _extractOrderedFaces(self, cell, faceOrderings):
        """Return faces for a 8-node hexahedron cell.
        """
        def orderingToFace(vertList):
            aFace = []
            for i in vertList:
                aFace.append(int(cell[i]))
            return aFace

        return [orderingToFace(o) for o in faceOrderings]

    def read(self):
        """
        0. Build cellsToVertices
        1. Recover needed vertexCoords and mapping from file using 
           cellsToVertices
        2. Build cellsToVertIDs proper from vertexCoords and vertex map
        3. Build faces
        4. Build cellsToFaces
        
        Isolate relevant data into three files, store in 
        `self.nodesFile` for $Nodes,
        `self.elemsFile` for $Elements. 
        `self.namesFile` for $PhysicalNames. 

        Returns vertexCoords, facesToVertexID, cellsToFaceID, 
                cellGlobalIDMap, ghostCellGlobalIDMap.
        """
        self.version, self.fileType, self.dataSize = self._getMetaData()
        self.nodesFile = self._isolateData("Nodes")
        self.elemsFile = self._isolateData("Elements")
        try:
            self.namesFile = self._isolateData("PhysicalNames")
        except EOFError, e:
            self.namesFile = None
            
        if self.dimensions is None:
            self.nodesFile.readline() # skip number of nodes

            # We assume we have a 2D file unless we find a node 
            # with a non-zero Z coordinate
            self.dimensions = 2
            for node in self.nodesFile:
                line   = node.split()
                
                newVert = [float(x) for x in line]
                
                if newVert[2] != 0.0:
                    self.dimensions = 3
                    break
                    
            self.nodesFile.seek(0)
            
        self.coordDimensions = self.coordDimensions or self.dimensions
            
        # we need a conditional here so we don't pick up 2D shapes in 3D
        if self.dimensions == 2: 
            self.numVertsPerFace = {1: 2, # 2-node line
                                    8: 2} # 3-node line
            self.numFacesPerCell = { 2: 3, # 3-node triangle (3 faces)
                                     9: 3, # 6-node triangle (we only read 1st 3)
                                    20: 3, # 9-node triangle (we only read 1st 3) 
                                    21: 3, # 10-node triangle (we only read 1st 3) 
                                    22: 3, # 12-node triangle (we only read 1st 3) 
                                    23: 3, # 15-node triangle (we only read 1st 3) 
                                    24: 3, # 15-node triangle (we only read 1st 3) 
                                    25: 3, # 21-node triangle (we only read 1st 3) 
                                     3: 4, # 4-node quadrangle (4 faces)
                                    10: 4, # 9-node quadrangle (we only read 1st 4)
                                    16: 4} # 8-node quadrangle (we only read 1st 4)
        elif self.dimensions == 3:
            self.numVertsPerFace = { 2: 3, # 3-node triangle (3 vertices)
                                     9: 3, # 6-node triangle (we only read 1st 3)
                                    20: 3, # 9-node triangle (we only read 1st 3) 
                                    21: 3, # 10-node triangle (we only read 1st 3) 
                                    22: 3, # 12-node triangle (we only read 1st 3) 
                                    23: 3, # 15-node triangle (we only read 1st 3) 
                                    24: 3, # 15-node triangle (we only read 1st 3) 
                                    25: 3, # 21-node triangle (we only read 1st 3) 
                                     3: 4, # 4-node quadrangle (4 vertices)
                                    10: 4, # 9-node quadrangle (we only read 1st 4)
                                    16: 4} # 8-node quadrangle (we only read 1st 4)
            self.numFacesPerCell = { 4: 4, # 4-node tetrahedron (4 faces)
                                    11: 4, # 10-node tetrahedron (we only read 1st 4)
                                    29: 4, # 20-node tetrahedron (we only read 1st 4)
                                    30: 4, # 35-node tetrahedron (we only read 1st 4)
                                    31: 4, # 56-node tetrahedron (we only read 1st 4)
                                     5: 6, # 8-node hexahedron (6 faces)
                                    12: 6, # 27-node tetrahedron (we only read 1st 6)
                                    17: 6, # 20-node tetrahedron (we only read 1st 6)
                                     6: 5, # 6-node prism (5 faces)
                                    13: 5, # 18-node prism (we only read 1st 6)
                                    18: 5, # 15-node prism (we only read 1st 6)
                                     7: 5, # 5-node pyramid (5 faces)
                                    14: 5, # 14-node pyramid (we only read 1st 5)
                                    19: 5} # 13-node pyramid (we only read 1st 5)
        else:
            raise GmshException("Mesh has fewer than 2 or more than 3 dimensions")

        print("Parsing elements.")
        (cellsData, 
         ghostsData, 
         facesData) = self._parseElementFile()
         
        self.cellsData = cellsData
        self.ghostsData = ghostsData
        self.facesData = facesData
        
        cellsToGmshVerts = cellsData.nodes + ghostsData.nodes
        numCellsTotal    = len(cellsData.nodes) + len(ghostsData.nodes)
        allShapeTypes    = cellsData.shapes + ghostsData.shapes
        allShapeTypes    = numpy.array(allShapeTypes)
        allShapeTypes    = numpy.delete(allShapeTypes, numpy.s_[numCellsTotal:])
        self.physicalCellMap = numpy.array(cellsData.physicalEntities 
                                        + ghostsData.physicalEntities)
        self.geometricalCellMap = numpy.array(cellsData.geometricalEntities 
                                           + ghostsData.geometricalEntities)

        if numCellsTotal < 1:
            errStr = "Gmsh hasn't produced any cells! Check your Gmsh code."
            errStr += "\n\nGmsh output:\n%s" % "".join(self.gmshOutput).rstrip()
            raise GmshException(errStr)

        print("Recovering coords.")
        print("numcells %d" % numCellsTotal)
        vertexCoords, vertIDtoIdx = self._vertexCoordsAndMap(cellsToGmshVerts)

        # translate Gmsh IDs to `vertexCoord` indices
        cellsToVertIDs = self._translateNodesToVertices(cellsToGmshVerts,
                                                        vertIDtoIdx)
        
        self.cellsToVerIds = cellsToVertIDs

        print("Building cells and faces.")
        (facesToV, 
         cellsToF,
         facesDict) = self._deriveCellsAndFaces(cellsToVertIDs,
                                                allShapeTypes,
                                                numCellsTotal)
            
        # cell entities were easy to record on parsing
        # but we don't use Gmsh faces, so we need to correlate the nodes 
        # that make up the Gmsh faces with the vertex IDs of the FiPy faces
        # so that we can check if any are named
        faceEntitiesDict = dict()
        
        # translate Gmsh IDs to `vertexCoord` indices
        facesToVertIDs = self._translateNodesToVertices(facesData.nodes,
                                                        vertIDtoIdx)
                                                    
        for face, physicalEntity, geometricalEntity in zip(facesToVertIDs, 
                                                           facesData.physicalEntities, 
                                                           facesData.geometricalEntities):
            faceEntitiesDict[' '.join([str(x) for x in sorted(face)])] = (physicalEntity, geometricalEntity)
            
        self.physicalFaceMap = numpy.zeros(facesToV.shape[-1:], 'l')
        self.geometricalFaceMap = numpy.zeros(facesToV.shape[-1:], 'l')
        for face in facesDict.keys():
            # not all faces are necessarily tagged
            if faceEntitiesDict.has_key(face):
                self.physicalFaceMap[facesDict[face]] = faceEntitiesDict[face][0]
                self.geometricalFaceMap[facesDict[face]] = faceEntitiesDict[face][1]
                
        self.physicalNames = self._parseNamesFile()
                                          
        print("Done with cells and faces.")
        return (vertexCoords, facesToV, cellsToF,  cellsToVertIDs,
                cellsData.idmap, ghostsData.idmap)
                
    def _vertexCoordsAndMap(self, cellsToGmshVerts):
        """
        Returns `vertexCoords` and mapping from Gmsh ID to `vertexCoords`
        indices (same as in MSHFile). 
        
        Unlike parent, doesn't use genfromtxt
        because we want to avoid loading the entire msh file into memory.
        """
        allVerts     = [v for c in cellsToGmshVerts for v in c] # flatten
        allVerts     = numpy.unique(numpy.array(allVerts, dtype=int)) # remove dups
        allVerts     = numpy.sort(allVerts)
        maxVertIdx   = allVerts[-1] + 1 # add one to offset zero
        vertGIDtoIdx = numpy.ones(maxVertIdx, 'l') * -1 # gmsh ID -> vertexCoords idx
        vertexCoords = numpy.empty((len(allVerts), self.coordDimensions))
        nodeCount    = 0

        # establish map. This works because allVerts is a sorted set.
        vertGIDtoIdx[allVerts] = numpy.arange(len(allVerts))

        self.nodesFile.readline() # skip number of nodes

        # now we walk through node file with a sorted unique list of vertices
        # in hand. When we encounter 0th element in `allVerts`, save it 
        # to `vertexCoords` then pop its ID off `allVerts`.
        for node in self.nodesFile:
            line   = node.split()
            nodeID = int(line[0])

            if nodeID == allVerts[nodeCount]:
                newVert = [float(x) for x in line[1:self.coordDimensions+1]]
                vertexCoords[nodeCount,:] = numpy.array(newVert)
                nodeCount += 1

            if len(allVerts) == nodeCount: 
                break

        # transpose for FiPy
        transCoords = vertexCoords.swapaxes(0,1)
        return transCoords, vertGIDtoIdx

    def _parseElementFile(self):
        """
        Return three objects, the first for non-ghost cells, the second for
        ghost cells, and the third for faces.

        All nastiness concerning ghost cell
        calculation is consolidated here: if we were ever to need to CALCULATE
        GHOST CELLS OURSELVES, the only code we'd have to change is in here.
        """
        cellsData = _ElementData()
        ghostsData = _ElementData()
        facesData = _ElementData()
        
        cellOffset = -1 # this will be subtracted from gmsh ID to
                        # obtain global ID
        faceOffset = -1 # this will be subtracted from gmsh ID to
                        # obtain global ID
                        
        pid = 1
        
        self.elemsFile.readline() # skip number of elements
        for el in self.elemsFile:
            currLineInts = [int(x) for x in el.split()]
            elemType     = currLineInts[1]

            if elemType in self.numFacesPerCell.keys():
                # element is a cell
                if cellOffset == -1:
                    # if first valid shape
                    cellOffset = currLineInts[0]
                currLineInts[0] -= cellOffset

                numTags = currLineInts[2]
                tags    = currLineInts[3:(3+numTags)]
                
                physicalEntity = tags.pop(0)
                geometricalEntity = tags.pop(0)

                if len(tags) > 0:
                    # next item is a count
                    if tags[0] != len(tags) - 1:
                        warnings.warn("Partition count %d does not agree with number of remaining tags %d." % (tags[0], len(tags) - 1), 
                                      SyntaxWarning, stacklevel=2)
                    tags = tags[1:]
            
                cellsData.add(currLine=currLineInts, elType=elemType, 
                              physicalEntity=physicalEntity, 
                              geometricalEntity=geometricalEntity)
            elif elemType in self.numVertsPerFace.keys():
                # element is a face
                if faceOffset == -1:
                    # if first valid shape
                    faceOffset = currLineInts[0]
                currLineInts[0] -= faceOffset

                numTags = currLineInts[2]
                tags    = currLineInts[3:(3+numTags)]
                
                # the partition tags for faces don't seem to always be present 
                # and don't always make much sense when they are
                
                physicalEntity = tags.pop(0)
                geometricalEntity = tags.pop(0)

                facesData.add(currLine=currLineInts, elType=elemType, 
                              physicalEntity=physicalEntity, 
                              geometricalEntity=geometricalEntity)
                              
        self.elemsFile.close() # tempfile trashed

        return cellsData, ghostsData, facesData

    def _parseNamesFile(self):
        physicalNames = {
            0: dict(),
            1: dict(),
            2: dict(),
            3: dict()
        }
        if self.namesFile is not None:
            self.namesFile.readline() # skip number of elements
            for nm in self.namesFile:
                nm = nm.split()
                if self.version > 2.0:
                    dim = [int(nm.pop(0))]
                else:
                    # Gmsh format prior to 2.1 did not unambiguously tie 
                    # physical names to physical entities of different dimensions
                    # http://article.gmane.org/gmane.comp.cad.gmsh.general/1601
                    dim = [0, 1, 2, 3]
                num = int(nm.pop(0))
                name = " ".join(nm)[1:-1]
                for d in dim:
                    physicalNames[d][name] = int(num)
                
        return physicalNames
        
    def makeMapVariables(self, mesh):
        """Utility function to make MeshVariables that define different domains in the mesh
        """
        from fipy.variables.cellVariable import CellVariable
        from fipy.variables.faceVariable import FaceVariable

        self.physicalCellMap = CellVariable(mesh=mesh, value=self.physicalCellMap)
        self.geometricalCellMap = CellVariable(mesh=mesh, value=self.geometricalCellMap)
        self.physicalFaceMap = FaceVariable(mesh=mesh, value=self.physicalFaceMap)
        self.geometricalFaceMap = FaceVariable(mesh=mesh, value=self.geometricalFaceMap)

        physicalCells = dict()
        for name in self.physicalNames[self.dimensions].keys():
            physicalCells[name] = (self.physicalCellMap == self.physicalNames[self.dimensions][name])
            
        physicalFaces = dict()
        for name in self.physicalNames[self.dimensions-1].keys():
            physicalFaces[name] = (self.physicalFaceMap == self.physicalNames[self.dimensions-1][name])
            
        return (self.physicalCellMap,
                self.geometricalCellMap,
                physicalCells,
                self.physicalFaceMap,
                self.geometricalFaceMap,
                physicalFaces)
                
class Gmsh2D():
    """Construct a 2D Mesh using Gmsh
    
    >>> radius = 5.
    >>> side = 4.
    >>> squaredCircle = Gmsh2D('''
    ... // A mesh consisting of a square inside a circle inside a circle
    ...                        
    ... // define the basic dimensions of the mesh
    ...                        
    ... cellSize = 1;
    ... radius = %(radius)g;
    ... side = %(side)g;
    ...                        
    ... // define the compass points of the inner circle
    ...                        
    ... Point(1) = {0, 0, 0, cellSize};
    ... Point(2) = {-radius, 0, 0, cellSize};
    ... Point(3) = {0, radius, 0, cellSize};
    ... Point(4) = {radius, 0, 0, cellSize};
    ... Point(5) = {0, -radius, 0, cellSize};
    ...                        
    ... // define the compass points of the outer circle
    ... 
    ... Point(6) = {-2*radius, 0, 0, cellSize};
    ... Point(7) = {0, 2*radius, 0, cellSize};
    ... Point(8) = {2*radius, 0, 0, cellSize};
    ... Point(9) = {0, -2*radius, 0, cellSize};
    ... 
    ... // define the corners of the square
    ... 
    ... Point(10) = {side/2, side/2, 0, cellSize/2};
    ... Point(11) = {-side/2, side/2, 0, cellSize/2};
    ... Point(12) = {-side/2, -side/2, 0, cellSize/2};
    ... Point(13) = {side/2, -side/2, 0, cellSize/2};
    ... 
    ... // define the inner circle
    ... 
    ... Circle(1) = {2, 1, 3};
    ... Circle(2) = {3, 1, 4};
    ... Circle(3) = {4, 1, 5};
    ... Circle(4) = {5, 1, 2};
    ... 
    ... // define the outer circle
    ... 
    ... Circle(5) = {6, 1, 7};
    ... Circle(6) = {7, 1, 8};
    ... Circle(7) = {8, 1, 9};
    ... Circle(8) = {9, 1, 6};
    ... 
    ... // define the square
    ... 
    ... Line(9) = {10, 13};
    ... Line(10) = {13, 12};
    ... Line(11) = {12, 11};
    ... Line(12) = {11, 10};
    ... 
    ... // define the three boundaries
    ... 
    ... Line Loop(1) = {1, 2, 3, 4};
    ... Line Loop(2) = {5, 6, 7, 8};
    ... Line Loop(3) = {9, 10, 11, 12};
    ... 
    ... // define the three domains
    ... 
    ... Plane Surface(1) = {2, 1};
    ... Plane Surface(2) = {1, 3};
    ... Plane Surface(3) = {3};
    ... 
    ... // label the three domains
    ... 
    ... // attention: if you use any "Physical" labels, you *must* label 
    ... // all elements that correspond to FiPy Cells (Physical Surace in 2D 
    ... // and Physical Volume in 3D) or Gmsh will not include them and FiPy
    ... // will not be able to include them in the Mesh. 
    ... 
    ... // note: if you do not use any labels, all Cells will be included.
    ... 
    ... Physical Surface("Outer") = {1};
    ... Physical Surface("Middle") = {2};
    ... Physical Surface("Inner") = {3};
    ... 
    ... // label the "north-west" part of the exterior boundary
    ... 
    ... // note: you only need to label the Face elements 
    ... // (Physical Line in 2D and Physical Surface in 3D) that correspond
    ... // to boundaries you are interested in. FiPy does not need them to
    ... // construct the Mesh.
    ... 
    ... Physical Line("NW") = {5};
    ... ''' % locals())

    It can be easier to specify certain domains and boundaries within Gmsh 
    than it is to define the same domains and boundaries with FiPy expressions.
    
    Here we compare obtaining the same Cells and Faces using FiPy's 
    parametric descriptions and Gmsh's labels.
    
    >>> x, y = squaredCircle.cellCenters

    >>> middle = ((x**2 + y**2 <= radius**2) 
    ...           & ~((x > -side/2) & (x < side/2)
    ...               & (y > -side/2) & (y < side/2)))

    >>> print (middle == squaredCircle.physicalCells["Middle"]).all()
    True
    
    >>> X, Y = squaredCircle.faceCenters

    >>> NW = ((X**2 + Y**2 > (1.99*radius)**2) 
    ...       & (X**2 + Y**2 < (2.01*radius)**2)
    ...       & (X <= 0) & (Y >= 0))

    >>> print (NW == squaredCircle.physicalFaces["NW"]).all()
    True
    """
    
    def __init__(self, filename, coordDimensions=2, order=1, background=None):
                     
        self.mshfile = MSHFile(filename=filename, 
                               dimensions=coordDimensions, 
                               coordDimensions=coordDimensions)
        
        (verts,
         faces,
         cells,
         _cells,
         self._cell_global_ids, 
         self._g_ell_global_ids) = self.mshfile.read()
         
        self.elements = numpy.array(_cells)

        self.mshfile.close()

        self.nfaces = faces.shape[-1]
        self.ncells = cells.shape[-1]        

        self.vertex_coords = verts
        self.face_vertex_ids = faces
        self.cell_face_ids = cells

        self._face_cell_ids = self._calcFaceCellIDs()

        # Transpose arrays
        self.face_cell_ids = self._face_cell_ids.filled(-1).T
        self.cell_face_ids = cells.T
        self.face_vertex_ids = faces.T
        
        # points is a NX3 array
        self.nnodes = self.vertex_coords.shape[-1]
        self.points = numpy.zeros( shape=(self.nnodes, 3) )
        self.points[:, 0] = self.vertex_coords[0]
        self.points[:, 1] = self.vertex_coords[1]

        self.type = MeshType.Gmsh

        # setup
        self._setup()

        print("Exiting Gmsh2D")

    def _setup(self):
        face_cell_ids = self.face_cell_ids
        cell_face_ids = self.cell_face_ids

        # boundary faces
        self.bfaces = bfaces = find(face_cell_ids == -1)/2

        # boundary cells
        self.bcells = bcells = {}
        for face in bfaces:
            cell = face_cell_ids[face][0]
            if not bcells.has_key(cell):
                bcells[cell] = [face]
            else:
                bcells[cell].append(face)

        # create ghost cells
        self._create_ghost_cells()

    def _create_ghost_cells(self):
        """Create ghost cells for the mesh"""
        elements = list(self.elements)
        points = list(self.points)
        face_vertex_ids = self.face_vertex_ids
        cell_face_ids = self.cell_face_ids
        face_cell_ids = self.face_cell_ids

        # get the boundary cells
        bcells = self.bcells

        nelements = self.ncells
        npoints = self.nnodes

        nghost = 0
        for cell in bcells:
            faces = bcells[cell] # face adjoining the boundary

            # vertices for the element 
            cell_nodes = elements[cell]

            for face in faces:
                face_nodes = face_vertex_ids[face]

                va = Vertex(*points[face_nodes[0]])
                vb = Vertex(*points[face_nodes[1]])

                tangent = Vector(va.x-vb.x,
                                 va.y-vb.y,
                                 va.z-vb.z)
                tangent.normalize()

                image_nodes = []
                for node in cell_nodes:
                    if node in face_nodes:
                        image_nodes.append(node)
                    else:
                        # now create the image point
                        p = Vertex(*points[node])
                        image = p.reflect(va, tangent)

                        # add the image point to the list of points
                        points.append( numpy.array([image.x, image.y, image.z]) )
                        image_nodes.append(npoints)
                        npoints += 1
                    
                # add the element nodes to the list of elements
                elements.append( image_nodes )

                # the ghost cell becomes the "right" cell
                face_cell_ids[face][1] = nelements
                nelements += 1
                nghost += 1

        self.nghost = nghost
        self.nelements = nelements
        self.npoints = npoints

        self.points = numpy.array( points )
        self.elements = numpy.array( elements )
        #self.elements = elements

    def _calcFaceCellIDs(self):
        array = MA.array(MA.indices(self.cell_face_ids.shape, 'l')[1], 
                         mask=MA.getmask(self.cell_face_ids))
        faceCellIDs = MA.zeros((2, self.nfaces), 'l')

        ## Nasty bug: MA.put(arr, ids, values) fills its ids and
        ## values arguments when masked!  This was not the behavior
        ## that was assumed when used below.  It was only working
        ## because the old fill value was 0 and the first element of
        ## the array needed to be 0 since the cell's face was
        ## 0. numerix.put() has been changed to deal with this
        ## properly.

        ##         MA.put(firstRow, cellFaceIDsFlat[::-1], array[::-1])
        ##         MA.put(secondRow, cellFaceIDsFlat, array)
        firstRow = faceCellIDs[0]
        secondRow = faceCellIDs[1]

        numpy.put(firstRow, self.cell_face_ids[::-1,::-1], array[::-1,::-1])
        numpy.put(secondRow, self.cell_face_ids, array)
        
        mask = ((False,) * self.nfaces, (firstRow == secondRow))
        return MA.sort(MA.array(faceCellIDs, mask = mask),
                                   axis=0)        

class _ElementData(object):
    """
    Bookkeeping for cells. Declared as own class for generality.

    "nodes": A Python list of the vertices that make up this element
    "shapes": A shapeTypes Python list
    "idmap": A Python list which maps vertexCoords idx -> global ID
    "physicalEntities": A Python list of the Gmsh physical entities each element is in
    "geometricalEntities": A Python list of the Gmsh geometrical entities each element is in
    """
    def __init__(self):
        self.nodes = []
        self.shapes = []
        self.idmap = [] # vertexCoords idx -> gmsh ID (global ID)
        self.physicalEntities = []
        self.geometricalEntities = []
        
    def add(self, currLine, elType, physicalEntity, geometricalEntity):
        numTags = currLine[2]
        self.nodes.append(currLine[(numTags+3):])
        self.shapes.append(elType)
        self.idmap.append(currLine[0])
        self.physicalEntities.append(physicalEntity)
        self.geometricalEntities.append(geometricalEntity)
